//*********************************************************************
//**                      3D geeomtry tools                          **
//**       Jarek Rossignac, October 2010, updates Oct 2011           **   
//**                 (points, vectors, and more)                     **   
//*********************************************************************

// ===== vector class
class vec { 
  float x=0, y=0, z=0; 
  vec () {
  }; 
  vec (float px, float py, float pz) {
    x = px; 
    y = py; 
    z = pz;
  };
  vec set (float px, float py, float pz) {
    x = px; 
    y = py; 
    z = pz; 
    return this;
  }; 
  vec set (vec V) {
    x = V.x; 
    y = V.y; 
    z = V.z; 
    return this;
  }; 
  vec add(vec V) {
    x+=V.x; 
    y+=V.y; 
    z+=V.z; 
    return this;
  };
  vec add(float s, vec V) {
    x+=s*V.x; 
    y+=s*V.y; 
    z+=s*V.z; 
    return this;
  };
  vec sub(vec V) {
    x-=V.x; 
    y-=V.y; 
    z-=V.z; 
    return this;
  };
  vec mul(float f) {
    x*=f; 
    y*=f; 
    z*=f; 
    return this;
  };
  vec div(float f) {
    x/=f; 
    y/=f; 
    z/=f; 
    return this;
  };
  vec div(int f) {
    x/=f; 
    y/=f; 
    z/=f; 
    return this;
  };
  vec rev() {
    x=-x; 
    y=-y; 
    z=-z; 
    return this;
  };
  float norm() {
    return(sqrt(sq(x)+sq(y)+sq(z)));
  }; 
  vec normalize() {
    float n=norm(); 
    if (n>0.000001) {
      div(n);
    }; 
    return this;
  };
  vec rotate(float a, vec I, vec J) {
    float x=d(this, I), y=d(this, J); 
    float c=cos(a), s=sin(a); 
    add(x*c-x-y*s, I); 
    add(x*s+y*c-y, J); 
    return this;
  }; // Rotate by a in plane (I,J)
} 
;

// ===== vector functions
vec V() {
  return new vec();
};                                                                          // make vector (x,y,z)
vec V(float x, float y, float z) {
  return new vec(x, y, z);
};                                            // make vector (x,y,z)
vec V(vec V) {
  return new vec(V.x, V.y, V.z);
};                                                          // make copy of vector V
vec A(vec A, vec B) {
  return new vec(A.x+B.x, A.y+B.y, A.z+B.z);
};                                       // A+B
vec A(vec U, float s, vec V) {
  return V(U.x+s*V.x, U.y+s*V.y, U.z+s*V.z);
};                               // U+sV
vec M(vec U, vec V) {
  return V(U.x-V.x, U.y-V.y, U.z-V.z);
};                                              // U-V
vec M(vec V) {
  return V(-V.x, -V.y, -V.z);
};                                                              // -V
vec V(vec A, vec B) {
  return new vec((A.x+B.x)/2.0, (A.y+B.y)/2.0, (A.z+B.z)/2.0);
}                      // (A+B)/2
vec V(vec A, float s, vec B) {
  return new vec(A.x+s*(B.x-A.x), A.y+s*(B.y-A.y), A.z+s*(B.z-A.z));
};      // (1-s)A+sB
vec V(vec A, vec B, vec C) {
  return new vec((A.x+B.x+C.x)/3.0, (A.y+B.y+C.y)/3.0, (A.z+B.z+C.z)/3.0);
};  // (A+B+C)/3
vec V(vec A, vec B, vec C, vec D) {
  return V(V(A, B), V(C, D));
};                                         // (A+B+C+D)/4
vec V(float s, vec A) {
  return new vec(s*A.x, s*A.y, s*A.z);
};                                           // sA
vec V(float a, vec A, float b, vec B) {
  return A(V(a, A), V(b, B));
}                                       // aA+bB 
vec V(float a, vec A, float b, vec B, float c, vec C) {
  return A(V(a, A, b, B), V(c, C));
}                   // aA+bB+cC
vec V(pt P, pt Q) {
  return new vec(Q.x-P.x, Q.y-P.y, Q.z-P.z);
};                                          // PQ
vec U(vec V) {
  float n = V.norm(); 
  if (n<0.000001) return V(0, 0, 0); 
  else return V(1./n, V);
};            // V/||V||
vec U(pt A, pt B) {
  return U(V(A, B));
}
vec N(vec U, vec V) {
  return V( U.y*V.z-U.z*V.y, U.z*V.x-U.x*V.z, U.x*V.y-U.y*V.x);
};                  // UxV CROSS PRODUCT (normal to both)
vec N(pt A, pt B, pt C) {
  return N(V(A, B), V(A, C));
};                                                   // normal to triangle (A,B,C), not normalized (proportional to area)
vec B(vec U, vec V) {
  return U(N(N(U, V), U));
}                                                           // (UxV)xV unit normal to U in the plane UV
vec R(vec V) {
  return V(-V.y, V.x, V.z);
} // rotated 90 degrees in XY plane
vec R(vec V, float a, vec I, vec J) {
  float x=d(V, I), y=d(V, J); 
  float c=cos(a), s=sin(a); 
  return A(V, V(x*c-x-y*s, I, x*s+y*c-y, J));
}; // Rotated V by a parallel to plane (I,J)

vec ratioVec(pt A, pt B, pt P) {
  vec va = U(V(A,B));
  vec vb = V(A,P);
  vec projToAB = V(dot(vb,va),va);
  return projToAB;
}

// ===== point class
class pt { 
  float x=0, y=0, z=0; 
  pt () {
  }; 
  pt (float px, float py, float pz) {
    x = px; 
    y = py; 
    z = pz;
  };
  pt set (float px, float py, float pz) {
    x = px; 
    y = py; 
    z = pz; 
    return this;
  }; 
  pt set (pt P) {
    x = P.x; 
    y = P.y; 
    z = P.z; 
    return this;
  }; 
  pt add(pt P) {
    x+=P.x; 
    y+=P.y; 
    z+=P.z; 
    return this;
  };
  pt add(vec V) {
    x+=V.x; 
    y+=V.y; 
    z+=V.z; 
    return this;
  };
  pt add(float s, vec V) {
    x+=s*V.x; 
    y+=s*V.y; 
    z+=s*V.z; 
    return this;
  };
  pt add(float dx, float dy, float dz) {
    x+=dx; 
    y+=dy; 
    z+=dz; 
    return this;
  };
  pt sub(pt P) {
    x-=P.x; 
    y-=P.y; 
    z-=P.z; 
    return this;
  };
  pt mul(float f) {
    x*=f; 
    y*=f; 
    z*=f; 
    return this;
  };
  pt mul(float dx, float dy, float dz) {
    x*=dx; 
    y*=dy; 
    z*=dz; 
    return this;
  };
  pt div(float f) {
    x/=f; 
    y/=f; 
    z/=f; 
    return this;
  };
  pt div(int f) {
    x/=f; 
    y/=f; 
    z/=f; 
    return this;
  };
  pt snap(float r) {
    float f=r/(sqrt(sq(x)+sq(y)+sq(z))); 
    x*=f; 
    y*=f; 
    z*=f; 
    return this;
  };
  pt coord2d(){
    return P(screenX(this.x, this.y, this.z),screenY(this.x, this.y, this.z),0); 
  }
  //   void projectOnCylinder(pt A, pt B, float r) {pt H = S(A,d(V(A,B),V(A,this))/d(V(A,B),V(A,B)),B); this.setTo(T(H,r,this));}
}
//  void projectOnCylinder(pt A, pt B, float r) {pt H = S(A,d(V(A,B),V(A,this))/d(V(A,B),V(A,B)),B); this.setTo(T(H,r,this));}   
// =====  point functions
pt P() {
  return new pt();
};                                            // point (x,y,z)
pt P(float x, float y, float z) {
  return new pt(x, y, z);
};                                            // point (x,y,z)
pt P(pt A) {
  return new pt(A.x, A.y, A.z);
};                                                           // copy of point P
pt P(pt A, float s, pt B) {
  return new pt(A.x+s*(B.x-A.x), A.y+s*(B.y-A.y), A.z+s*(B.z-A.z));
};        // A+sAB
pt P(pt A, pt B) {
  return P((A.x+B.x)/2.0, (A.y+B.y)/2.0, (A.z+B.z)/2.0);
}                             // (A+B)/2
pt P(pt A, pt B, pt C) {
  return new pt((A.x+B.x+C.x)/3.0, (A.y+B.y+C.y)/3.0, (A.z+B.z+C.z)/3.0);
};     // (A+B+C)/3
pt P(pt A, pt B, pt C, pt D) {
  return P(P(A, B), P(C, D));
};                                            // (A+B+C+D)/4
pt P(float s, pt A) {
  return new pt(s*A.x, s*A.y, s*A.z);
};                                            // sA
pt A(pt A, pt B) {
  return new pt(A.x+B.x, A.y+B.y, A.z+B.z);
};                                         // A+B
pt P(float a, pt A, float b, pt B) {
  return A(P(a, A), P(b, B));
}                                        // aA+bB 
pt P(float a, pt A, float b, pt B, float c, pt C) {
  return A(P(a, A), P(b, B, c, C));
}                     // aA+bB+cC 
pt P(float a, pt A, float b, pt B, float c, pt C, float d, pt D) {
  return A(P(a, A, b, B), P(c, C, d, D));
}   // aA+bB+cC+dD
pt P(pt P, vec V) {
  return new pt(P.x + V.x, P.y + V.y, P.z + V.z);
}                                 // P+V
pt P(pt P, float s, vec V) {
  return new pt(P.x+s*V.x, P.y+s*V.y, P.z+s*V.z);
}                           // P+sV
pt P(pt O, float x, vec I, float y, vec J) {
  return P(O.x+x*I.x+y*J.x, O.y+x*I.y+y*J.y, O.z+x*I.z+y*J.z);
}  // O+xI+yJ
pt P(pt O, float x, vec I, float y, vec J, float z, vec K) {
  return P(O.x+x*I.x+y*J.x+z*K.x, O.y+x*I.y+y*J.y+z*K.y, O.z+x*I.z+y*J.z+z*K.z);
}  // O+xI+yJ+kZ
pt R(pt P, float a, vec I, vec J, pt G) {
  float x=d(V(G, P), I), y=d(V(G, P), J); 
  float c=cos(a), s=sin(a); 
  return P(P, x*c-x-y*s, I, x*s+y*c-y, J);
}; // Rotated P by a around G in plane (I,J)
void makePts(pt[] C) {
  for (int i=0; i<C.length; i++) C[i]=P();
} // fills array C with points initialized to (0,0,0)
pt Predict(pt A, pt B, pt C) {
  return P(B, V(A, C));
};     // B+AC, parallelogram predictor
void v(pt P) {
  vertex(P.x, P.y, P.z);
} // rendering


// ===== mouse tools
pt Mouse() {
  return P(mouseX, mouseY, 0);
};                                          // current mouse location
pt Pmouse() {
  return P(pmouseX, pmouseY, 0);
};
vec MouseDrag() {
  return V(mouseX-pmouseX, mouseY-pmouseY, 0);
};                     // vector representing recent mouse displacement

// ===== measures
float d(vec U, vec V) {
  return U.x*V.x+U.y*V.y+U.z*V.z;
};                                            //U*V dot product
float dot(vec U, vec V) {
  return d(U,V);
};
float d(pt P, pt Q) {
  return sqrt(sq(Q.x-P.x)+sq(Q.y-P.y)+sq(Q.z-P.z));
};                            // ||AB|| distance
float d2(pt P, pt Q) {
  return sq(Q.x-P.x)+sq(Q.y-P.y)+sq(Q.z-P.z);
};                                 // AB^2 distance squared
float m(vec U, vec V, vec W) {
  return d(U, N(V, W));
};                                                 // (UxV)*W  mixed product, determinant
float m(pt E, pt A, pt B, pt C) {
  return m(V(E, A), V(E, B), V(E, C));
}                                    // det (EA EB EC) is >0 when E sees (A,B,C) clockwise
float n2(vec V) {
  return sq(V.x)+sq(V.y)+sq(V.z);
};                                                   // V*V    norm squared
float n(vec V) {
  return sqrt(n2(V));
};                                                                // ||V||  norm
float area(pt A, pt B, pt C) {
  return n(N(A, B, C))/2;
};                                               // area of triangle 
float volume(pt A, pt B, pt C, pt D) {
  return m(V(A, B), V(A, C), V(A, D))/6;
};                           // volume of tet 
boolean parallel (vec U, vec V) {
  return n(N(U, V))<n(U)*n(V)*0.00001;
}                               // true if U and V are almost parallel
float angle(vec U, vec V) {
  return acos(d(U, V)/n(V)/n(U));
};                                         // angle(U,V)
boolean cw(vec U, vec V, vec W) {
  return m(U, V, W)>=0;
};                                              // (UxV)*W>0  U,V,W are clockwise
boolean cw(pt A, pt B, pt C, pt D) {
  return volume(A, B, C, D)>=0;
};                                    // tet is oriented so that A sees B, C, D clockwise 

// ===== rotate 

// ===== render
void normal(vec V) {
  normal(V.x, V.y, V.z);
};                                          // changes normal for smooth shading
void vertex(pt P) {
  vertex(P.x, P.y, P.z);
};                                           // vertex for shading or drawing
void vTextured(pt P, float u, float v) {
  vertex(P.x, P.y, P.z, u, v);
};                          // vertex with texture coordinates
void show(pt P, pt Q) {
  line(Q.x, Q.y, Q.z, P.x, P.y, P.z);
};                       // draws edge (P,Q)
void show(pt P, vec V) {
  line(P.x, P.y, P.z, P.x+V.x, P.y+V.y, P.z+V.z);
};          // shows edge from P to P+V
void show(pt P, float d, vec V) {
  line(P.x, P.y, P.z, P.x+d*V.x, P.y+d*V.y, P.z+d*V.z);
}; // shows edge from P to P+dV
void show(pt A, pt B, pt C) {
  beginShape(); 
  vertex(A);
  vertex(B); 
  vertex(C); 
  endShape(CLOSE);
};                      // volume of tet 
void show(pt A, pt B, pt C, pt D) {
  beginShape(); 
  vertex(A); 
  vertex(B); 
  vertex(C); 
  vertex(D); 
  endShape(CLOSE);
};                      // volume of tet 
void show(pt P, float r) {
  pushMatrix(); 
  translate(P.x, P.y, P.z); 
  sphere(r); 
  popMatrix();
}; // render sphere of radius r and center P

void show(pt P, float s, vec I, vec J, vec K) {
  noStroke(); 
  fill(yellow); 
  show(P, 5); 
  stroke(red); 
  show(P, s, I); 
  stroke(green); 
  show(P, s, J); 
  stroke(blue); 
  show(P, s, K);
}; // render sphere of radius r and center P
void show(pt P, String s) {
  text(s, P.x, P.y, P.z);
}; // prints string s in 3D at P
void show(pt P, String s, vec D) {
  text(s, P.x+D.x, P.y+D.y, P.z+D.z);
}; // prints string s in 3D at P+D

// ==== curve
void bezier(pt A, pt B, pt C, pt D) {
  bezier(A.x, A.y, A.z, B.x, B.y, B.z, C.x, C.y, C.z, D.x, D.y, D.z);
} // draws a cubic Bezier curve with control points A, B, C, D
void bezier(pt [] C) {
  bezier(C[0], C[1], C[2], C[3]);
} // draws a cubic Bezier curve with control points A, B, C, D
pt bezierPoint(pt[] C, float t) {
  return P(bezierPoint(C[0].x, C[1].x, C[2].x, C[3].x, t), bezierPoint(C[0].y, C[1].y, C[2].y, C[3].y, t), bezierPoint(C[0].z, C[1].z, C[2].z, C[3].z, t));
}
vec bezierTangent(pt[] C, float t) {
  return V(bezierTangent(C[0].x, C[1].x, C[2].x, C[3].x, t), bezierTangent(C[0].y, C[1].y, C[2].y, C[3].y, t), bezierTangent(C[0].z, C[1].z, C[2].z, C[3].z, t));
}
void PT(pt P0, vec T0, pt P1, vec T1) {
  float d=d(P0, P1)/3;  
  bezier(P0, P(P0, -d, U(T0)), P(P1, -d, U(T1)), P1);
} // draws cubic Bezier interpolating  (P0,T0) and  (P1,T1) 
void PTtoBezier(pt P0, vec T0, pt P1, vec T1, pt [] C) {
  float d=d(P0, P1)/3;  
  C[0].set(P0); 
  C[1].set(P(P0, -d, U(T0))); 
  C[2].set(P(P1, -d, U(T1))); 
  C[3].set(P1);
} // draws cubic Bezier interpolating  (P0,T0) and  (P1,T1) 
vec vecToCubic (pt A, pt B, pt C, pt D, pt E) {
  return V( (-A.x+4*B.x-6*C.x+4*D.x-E.x)/6, (-A.y+4*B.y-6*C.y+4*D.y-E.y)/6, (-A.z+4*B.z-6*C.z+4*D.z-E.z)/6);
}
vec vecToProp (pt B, pt C, pt D) {
  float cb=d(C, B);  
  float cd=d(C, D); 
  return V(C, P(B, cb/(cb+cd), D));
};  

// ==== perspective
pt Pers(pt P, float d) { 
  return P(d*P.x/(d+P.z), d*P.y/(d+P.z), d*P.z/(d+P.z) );
};

pt InverserPers(pt P, float d) { 
  return P(d*P.x/(d-P.z), d*P.y/(d-P.z), d*P.z/(d-P.z) );
};

// ==== intersection
boolean intersect(pt P, pt Q, pt A, pt B, pt C, pt X) {
  return intersect(P, V(P, Q), A, B, C, X);
} // if (P,Q) intersects (A,B,C), return true and set X to the intersection point

boolean intersect(pt E, vec T, pt A, pt B, pt C, pt X) { // if ray from E along T intersects triangle (A,B,C), return true and set X to the intersection point
  vec EA=V(E, A), EB=V(E, B), EC=V(E, C), AB=V(A, B), AC=V(A, C); 
  boolean s=cw(EA, EB, EC), sA=cw(T, EB, EC), sB=cw(EA, T, EC), sC=cw(EA, EB, T); 
  if ( (s==sA) && (s==sB) && (s==sC) ) return false;
  float t = m(EA, AC, AB) / m(T, AC, AB);
  X.set(P(E, t, T));
  return true;
}

boolean rayIntersectsTriangle(pt E, vec T, pt A, pt B, pt C) { // true if ray from E with direction T hits triangle (A,B,C)
  vec EA=V(E, A), EB=V(E, B), EC=V(E, C); 
  boolean s=cw(EA, EB, EC), sA=cw(T, EB, EC), sB=cw(EA, T, EC), sC=cw(EA, EB, T); 
  return  (s==sA) && (s==sB) && (s==sC) ;
}

boolean edgeIntersectsTriangle(pt P, pt Q, pt A, pt B, pt C) {
  vec PA=V(P, A), PQ=V(P, Q), PB=V(P, B), PC=V(P, C), QA=V(Q, A), QB=V(Q, B), QC=V(Q, C); 
  boolean p=cw(PA, PB, PC), q=cw(QA, QB, QC), a=cw(PQ, PB, PC), b=cw(PA, PQ, PC), c=cw(PQ, PB, PQ); 
  return (p!=q) && (p==a) && (p==b) && (p==c);
}

float rayParameterToIntersection(pt E, vec T, pt A, pt B, pt C) {
  vec AE=V(A, E), AB=V(A, B), AC=V(A, C); 
  return - m(AE, AC, AB) / m(T, AC, AB);
}

float angleDraggedAround(pt G) {  // returns angle in 2D dragged by the mouse around the screen projection of G
  pt S=P(screenX(G.x, G.y, G.z), screenY(G.x, G.y, G.z), 0);
  vec T=V(S, Pmouse()); 
  vec U=V(S, Mouse());
  return atan2(d(R(U), T), d(U, T));
}

float toRad(float a) {
  return(a*PI/180);
}      
int toDeg(float a) {
  return int(a*180/PI);
}  

void showShrunkOffset(pt A, pt B, pt C, float e, float h) {
  vec N=U(N(V(A, B), V(A, C))); 
  showShrunk(P(A, h, N), P(B, h, N), P(C, h, N), e);
} // offset by h along normal

void showShrunk(pt A, pt B, pt C, float e) { // shrink by e
  vec AB = U(V(A, B)), BC = U(V(B, C)), CA = U(V(C, A));
  float a = e/n(N(CA, AB)), b = e/n(N(AB, BC)), c = e/n(N(BC, CA));
  float d = max(d(A, B)/3, d(B, C)/3, d(C, A)/3);
  a=min(a, d); 
  b=min(b, d); 
  c=min(c, d);
  pt As=P(A, a, AB, -a, CA), Bs=P(B, b, BC, -b, AB), Cs=P(C, c, CA, -c, BC);
  beginShape(); 
  vertex(As); 
  vertex(Bs); 
  vertex(Cs); 
  endShape(CLOSE);
} 

float scaleDraggedFrom(pt G) {
  pt S=P(screenX(G.x, G.y, G.z), screenY(G.x, G.y, G.z), 0); 
  return d(S, Mouse())/d(S, Pmouse());
}

// INTERPOLATING CURVE
void drawCurve(pt A, pt B, pt C, pt D) {
  float d=d(A, B)+d(B, C)+d(C, D); 
  beginShape(); 
  for (float t=0; t<=1; t+=0.025) vertex(P(A, B, C, D, t*d)); 
  endShape();
}
void drawSamplesOnCurve(pt A, pt B, pt C, pt D, float r) {
  float d=d(A, B)+d(B, C)+d(C, D); 
  for (float t=0; t<=1; t+=0.025) show(P(A, B, C, D, t*d), r);
}
pt P(pt A, pt B, pt C, pt D, float t) {
  return P(0, A, d(A, B), B, d(A, B)+d(B, C), C, d(A, B)+d(B, C)+d(C, D), D, t);
}
pt P(float a, pt A, float b, pt B, float c, pt C, float d, pt D, float t) {
  pt E = P(A, (t-a)/(b-a), B), F = P(B, (t-b)/(c-b), C), G = P(C, (t-c)/(d-c), D), 
  H = P(E, (t-a)/(c-a), F), I = P(F, (t-b)/(d-b), G);
  return P(H, (t-a)/(d-a), I);
}

pt NUBS(float a, pt A, float b, pt B, float c, pt C, float d, pt D, float e, float t) {
  pt E = P(A, (a+b+t*c)/(a+b+c), B), F = P(B, (b+t*c)/(b+c+d), C), G = P(C, (t*c)/(c+d+e), D), 
  H = P(E, (b+t*c)/(b+c), F), I = P(F, (t*c)/(c+d), G), 
  J = P(H, t, I);
  return J;
}

// Linear
pt L(pt A, pt B, float t) {
  return P(A.x+t*(B.x-A.x), A.y+t*(B.y-A.y), A.z+t*(B.z-A.z));
}
pt L3(pt A, pt B, pt C, float t){ //t is 0-1 and 1-2
  return L(L(A,B,t),L(B,C,t-1),t/2);
}
pt L4(pt A, pt B, pt C, pt D, float t){ //t is 0-1 between B and C
  pt pOne = L3(A,B,C,1+t);
  pt pTwo = L3(B,C,D,t);
  return P(pOne, pTwo); //P(pt, pt) is the midpoint formula
}


// Interpolation non-uniform (Neville's algorithm)
pt I(float a, pt A, float b, pt B, float t) {
  return L(A, B, (t-a)/(b-a));
}                               // P(a)=A, P(b)=B
pt I(float a, pt A, float b, pt B, float c, pt C, float t) {
  pt P=I(a, A, b, B, t); 
  pt Q=I(b, B, c, C, t); 
  return I(a, P, c, Q, t);
} // P(a)=A, P(b)=B, P(c)=C
pt I(float a, pt A, float b, pt B, float c, pt C, float d, pt D, float t) {
  pt P=I(a, A, b, B, c, C, t); 
  pt Q=I(b, B, c, C, d, D, t); 
  return I(a, P, d, Q, t);
} // P(a)=A, P(b)=B, P(c)=C, P(d)=D
pt I(float a, pt A, float b, pt B, float c, pt C, float d, pt D, float e, pt E, float t) {
  pt P=I(a, A, b, B, c, C, d, D, t); 
  pt Q=I(b, B, c, C, d, D, e, E, t); 
  return I(a, P, e, Q, t);
}
// Interpolation proportional to Distance
pt D(pt A, pt B, pt C, float t) {
  float a=0, b=d(A, B), c=b+d(B, C); 
  return I(a, A, b, B, c, C, a+t*(c-a));
}
pt D(pt A, pt B, pt C, pt D, float t) {
  float a=0, b=d(A, B), c=b+d(B, C), d=c+d(C, D); 
  return I(a, A, b, B, c, C, d, D, a+t*(d-a));
}
pt D(pt A, pt B, pt C, pt D, pt E, float t) {
  float a=0, b=d(A, B), c=b+d(B, C), d=c+d(C, D), e=d+d(D, E); 
  return I(a, A, b, B, c, C, d, D, e, E, a+t*(e-a));
}

